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For the past decade there has been a considerable debate about the existence of chaos in the mixmaster cosmological model. 
The debate has been hampered by the coordinate, or observer dependence of standard chaotic indicators such as Lyapunov 
exponents. Here we use coordinate independent, fractal methods to show the mixmaster universe is indeed chaotic. 
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The origin of the universe and the fate of collapsing 
stars are two of the great mysteries in nature. In general 
relativity without exotic matter, the singularity theorems 
of Hawking and Penrose |^ argue that the gravitational 
collapse of very massive stars ends singular and that the 
universe was born singular. These singular settings force 
gravity to face quantum mechanics. As well as expos- 
ing the fundamental laws of physics, the singular cores 
of black holes and the origin of the cosmos draw deep 
connections to the laws of thermodynamics and, as we 
will discuss here, to chaos. 

Earlier, Khalatnikov and Lifshitz j2j argued that sin- 
gular solutions were the exception rather than the rule, 
putting them at odds with the singularity theorems. 
Their argument was that deformations in spacetime 
would be amplified during collapse and consequently 
would fight the formation of a singularity. This implied 
that the known symmetric singular solutions were atyp- 
ical. The conflict was resolved when they realized the 
singularity in a collapsed star could be chaotic [H . They 
conjectured that a generic singularity drives spacetime to 
churn and oscillate chaotically. Independently, Misner ^] 
suggested a chaotic approach to an early universe singu- 
larity. In his mixmaster universe, the different directions 
in 3-space alternate in cycles of anisotropic collapse and 
expansion. A popular account of these developments may 
be found in Thome's recent book Q. 

While the emergence of chaos helped our understand- 
ing of singularities in general relativity, building a re- 
silient theory of relativistic chaos has become a task of 
its own. A debate has raged over whether or not the 
mixmaster universe is chaotic. Studies of the mixmaster 
dynamics using both approximate maps and numer- 
ical integrations have each yielded conflicting results 
as to the existence of positive Lyapunov exponents - a 
standard chaotic indicator. Finally, it was realised PpO| 
that Lyapunov exponents are coordinate dependent and 
the conflicting results were a consequence of the differ- 
ent coordinate systems. In short, Lyapunov exponents 
are not reliable indicators of chaos in general relativity. 
Using a different approach, it was shown that the mix- 
master equations fail the Painleve test |jll[ . This suggests 
that the mixmaster may be chaotic, but the Painleve test 



is also inconclusive. A detailed review of the mixmaster 
debate can be found in Ref. fl^ . 

In this letter we show that the mixmaster universe is 
indeed chaotic by using coordinate independent, fractal 
methods. A fractal set of self-similar universes is uncov- 
ered by numerically solving Einstein's equations. These 
universes form fractal boundaries in the space of ini- 
tial conditions. Such fractal partitions are the result of 
chaotic dynamics. We emphasize that our approach can 
be used to study any system described by general relativ- 
ity. The mixmaster is studied here as a topical example. 

The mixmaster universe Q has closed spatial sections 
with the topology of a three sphere. The vacuum field 
equations of general relativity lead to the equations 

(In a^)" = (62 _ c2)2 _ a4 ^^^^ 5^ 

Here (a, 6, c) are the scale factors for the three spatial 
axes, a prime denotes d/dr and dt = abcdr. In a numer- 
ical study it is advantageous to use a = Ina, /3 — Inb, 
7 = Inc and dT = dt / {abcla{abc)) as integration vari- 
ables. These variables cautiously approach the curvature 
singularity ai t — 0. Asymptotically the new time vari- 
able is related to the comoving time thy T — ln(ln(l/i)). 
The equations of motion (|l|) can be integrated to yield 
the Hamiltonian constraint H = 0, where 

H ^ (lna)'(ln6)' + (Ina)'(lnc)' + (ln6)'(lnc)' 

-- [a* + 6^ + c4 - 2 {aH^ + b^c" + a^c^)] . 



(2) 



When the potential terms on the r.h.s. of (||) are small, 
the mixmaster coasts in an approximate Kasner phase 
described by 



ds^ 



-dt^ 



t^Padx^+t^Pbdy^ + t^Pcdz\ 
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The indices can be written as pi — (1 + u)/(l + u + w^), 
P2 — -~u/{l+u+u'^), andp3 = {u+u'^)/ {1+u+u'^), where 
(a, 5, c) can take any ordering of (1, 2, 3) and u £ [1, 00). 

Initial conditions can be set on the surface a — 0, 
da/dT = a < 0, with the four variables {u, v, s, S) p3|: 



a = , 
a = sT,pi 



/3 = 



1 + V + uv 
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I](w + uv) 
1 + V + uv 
sEp3 . 
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The anisotropy in the sizes and velocities of the axes 
is quantified by v and u. The variables S and s are 
overall scale factors. We will primarily be interested in 
the evolution of (m, u) as the mixmaster singularity is 
approached. That is, we are mostly interested in the 
relative rates of expansion of the three spatial directions. 

Before turning to the full mixmaster dynamics, we con- 
sider the properties of a two dimensional map that ap- 
proximately describes the evolution of {u, v) as the sin- 
gularity is approached. We use the map to guide our 
study of the full dynamics. A detailed comparison of our 
treatment and the original studies of the related Gauss 
map will be delineated elsewhere We men- 

tion that unlike previous treatments, we focus on chaotic 
scattering. A complete description of chaotic scattering 
is encoded in the unstable periodic orbits |]l5|,^. These 
orbits form what is know as a strange repellor - a close 
cousin of the more familiar strange attractor. We expose 
the fractal nature of the repellor. The same fractal set is 
then found in the full dynamics. 

The map evolves {u, v) forward in discrete jumps |p|,|l3|, 

u— 1, V + 1 , u > 2 (oscillations) 
F{u,v)^{ 1 1 (5) 



u — 1 ' V 



1 , u < 2 (bounces) . 



During an oscillation, one pair of axes oscillates out of 
phase while the third decreases monotonically. At a 
bounce, the roles of the three axes are interchanged and 
a different axis decreases monotonically. 

The strange repellor is the fractal set of points along 
periodic orbits, F''{u,v) = {u,v). Physically these or- 
bits are self-similar universes. Periodic orbits with pe- 
riod p < k can be divided into m oscillations and k — m 
bounces. Since bounces and oscillations do not commute, 
the number of u fixed points on an orbit with m oscil- 
lations is ''Cm- The total number of fixed points at or- 
der k is given by the sum over all possibilities squared, 
N{k) = (2'= - 1)2. Thus, the topological entropy of 
the strange repellor is given by 



lim Y\nN{k) 

fe^oo k 
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Since Ht > the i^-map is chaotic. 

To make contact with the continuum dynamics, we 
display a portion of the repellor's future invariant set 
in Fig. 1. The future invariant set is the collection of 
lines F'^(u) = u, {v arbitary). The sequence of gaps 
around the rationals form what is known as a Farey 
tree. In a complementary fashion, the repellor lies on 
the periodic irrationals of the irrational Farey tree [ p4[ . 
A similar collection of lines occurs in each integer interval 
u = [n,n + l], but with exponentially decreasing density. 
The fall off can also be understood from the combina- 
torics of oscillations O and bounces B. A root in this 



interval corresponds to the word, 0("~^^i3. The number 
of n-letter words that can be formed from a 2-letter al- 
phabet is 2", so the fraction of roots in each interval is 
2~". The map and the continuum dynamics will be seen 
to describe the same future invariant set. 

The universes which comprise the strange repellor for- 
ever repeat some prescribed cycle in (u^v). In contrast, 
a typical mixmaster universe will launch into an infinite 
pattern of oscillations and bounces that never repeats. 
Moreover, typical universes have an invariant probability 
~ u/ \nu that grows with u [ pT| . As a result, typical uni- 
verses scatter to m — > cxd while universes on the repellor 
are concentrated at small u values. 
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FIG. 1. The map's future invariant set in the interval 
1 < It < 2. 

We can quantify the multifractal nature of the repellor 
by measuring its fractal dimensions Dq, where 



lim.l^^-i' 



In € 



(7) 



Here N{£) are the number of hypercubes of side length e 
needed to cover the fractal and pi is the fraction of points 
in the i*"^ hypercube. The standard capacity dimension is 
recovered when q — 0, the information dimension when 
q = 1, etc. For homogeneous fractals all the various 
dimensions yield the same result. The dimensions Dq 
are invariant under diffeomorphisms for all q. 

Since the periodic orbits of (^ are everywhere dense, it 
follows that the future invariant set in Fig. 1 has Do = 2. 
However, points on a small period orbit are visited with 
greater frequency and so have a larger pi than high pe- 
riod orbits. This generates an uneven distribution which 
ensures that Dq < 2 ior q > 0. Numerically solving for 
all roots up to fc = 16 we find Di = 1.87 ± 0.01, thus 
confirming the multifractal nature of the repellor. 

If the F-map had been obtained from the full Einstein 
equations without approximation, we could conclude that 
the mixmaster universe is chaotic. Since approximations 
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were made [^,|3|, the possibilty remains that the full 
equations are integrable. The approximations may have 
failed to preserve some integrals of the motion, thus lead- 
ing to a false chaotic signal. We show this is not the case. 

Since the mixmaster phase space is not compact, any 
chaotic behaviour is likely to be transient. There is a 
standard procedure to investigate such chaotic scatter- 
ing. First we identify the different asymptotic outcomes 
the system might have. Once outcomes are identified, 
several methods can be used to search for a strange re- 
pellor. The simplest method looks for a fractal pattern in 
plots of an appropriately defined scattering angle and im- 
pact parameter. Alternatively, the strange repellor can 
be hunted directly with a numerical shooting procedure 
called PIM triples 

Our prefered method is to look for fractal basin bound- 
aries . Each outcome has a basin of attraction in the 
space of initial conditions. We may plot these basins by 
assigning a different color to each outcome, and then col- 
oring all initial conditions according to their outcome. If 
the boundaries between these outcomes are smooth, then 
the dynamics is regular. Conversely, if the boundaries 
are fractal, the dynamics is chaotic. The set of points 
belonging to the fractal boundary form the strange re- 
pellor's future invariant set. 

The power of these methods as a tool for studying 
chaos in general relativity is twofold. First, a fractal is 
a non-differentiable structure and so cannot be removed 
by any differentiable coordinate transformation. Thus, 
a fractal basin boundary provides an observer indepen- 
dent signal of chaos. Second, most systems in general 
relativity have non-compact phase spaces, so most chaos 
will be transient. Other coordinate independent meth- 
ods of studying chaos in general relativity, such as curva- 
ture based methods |21|, only work for compact systems. 
Previously fractal techniques were used to show there is 
chaos in multi- black hole spacetimes [ p2| , and in various 
inflationary cosmological models [2^ ]. 

As it stands, the mixmaster dynamics rarely leads to 
definite outcomes. For typical trajectories the sequence 
of oscillations and bounces continues ad infinitum as the 
singularity is approached. A trajectory that visits any 
finite value of u, no matter how large, will eventually re- 
turn to bounce again. The way around this asymptotic 
backwash problem is to assign outcomes at some large, 
but finite distance from the scattering region. This proce- 
dure cannot lead to a false chaotic signature. At worst, 
it will return smooth boundaries for a chaotic system. 
This occurs when trajectories near the strange repellor 
are prematurely assigned to a particular outcome. 

For the mixmaster universe, typical trajectories scat- 
ter to large u values, while orbits on the repellor are 
concentrated at small u values. Accordingly, we assign 
outcomes when u becomes larger than some fixed num- 
ber, Mmax, during a Kasner phase. Due to the system's 
SO{2>) symetry, a large u value leads to three equally 



likely outcomes. Physically, these are the three states 
of highly anisotropic expansion with pi w p2 ~ and 
P3 w 1. Thus, the space of initial conditions can be color 
coded depending on which axis is collapsing most quickly. 
We color these black for a, grey for h and white for c. 




FIG. 2. Basin boundaries in the (u, v) plane for the full, 
unapproximated mixmaster dynamics. 

Our prescription is very easy to implement numerically. 
In order to faciliate comparison with Fig. 1, we chose ini- 
tial condition in accordance with (^) by selecting (uo, Vq) 
from a 300^ grid, setting So = 1 and then using (||) to fix 
Sq. The initial conditions are then evolved according to 
the equations of motion (Q), and the ratios of a, /? and 
7 are monitored to see if the universe is in an approxi- 
mate Kasner phase described by (H). If it is, the value of 
u is extracted and an outcome is assigned if u > Umax- 
We choose Mmax = 7 since there is only a 1 in 2^ chance 
that a trajectory with u > 7 lies on the strange repel- 
lor. Moreover, typical aperiodic trajectories will scatter 
to u > Winax after a few bounces, so the numerical inte- 
gration is kept short and numerical errors do not become 
large. To confirm this, the unenforced Hamiltonian con- 
straint was monitored at all times and found to be 
satisfied within numerical tolerances. 

In Fig. 2 we display the basin boundaries in a portion 
of the (u, v) plane. We see a complicated mixture of both 
regular and fractal basin boundaries. The numerically 
generated basin boundaries are built of universes which 
ride the repellor for many orbits before being thrown off. 
Similar fractal basins can be found by viewing alternative 
slices through phase space, such as the {/3, $) plane. The 
overall morphology of the basins is altered little by de- 
manding more strongly anisotropic outcomes. From Fig. 
3 we see that the fractal nature of the boundary persists 
on finer and finer scales. 

Aside from some mild warpage, the future invariant 
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set (basin boundaries) seen in Figs. 2 & 3 for the full dy- 
namics is strikingly similar to that shown in Fig. 1 for the 
discrete map. The warpage can be accounted for by our 
choice of initial conditions near the maximum of expan- 
sion, where the approximations used to derive the map 
break down. However, the important fine scale structure 
is laid down much nearer the singularity, where the map 
works well, and this accounts for the agreement in the 
detailed structure seen in Figs. 1 and 3. Indeed, a calcu- 
lation of the information dimension of Fig. 3 using 
the uncertainty exponent method yields Di = 1.86±0.01, 
in agreement with the map. 



ities are chaotic. This in turn suggests that quantum 
gravity may have to confront quantum chaos. 

We thank J. Barrow, C. Dcttmann, N. Frankel and 
P. Ferreira for helpful comments. We are grateful to C. 
Dettmann for letting us adapt his computer codes. 




FIG. 3. A portion of Fig. 2 magnified 25 times. 

By exploiting techniques originally developed to study 
chaotic scattering, we gain a new perspective on the evo- 
lution of the mixmaster cosmology. We found a fractal 
structure, the strange repellor, describes the chaos well. 
The strange repellor is the collection of all universes pe- 
riodic in {u,v). A typical, aperiodic universe will expe- 
rience a transient age of chaos if it brushes against the 
repellor. The fractal was exposed in both the exact Ein- 
stein equations and the discrete map used to approximate 
the evolution. Most importantly, our fractal approach is 
independent of the time coordinate used. An outcome is 
an outcome no matter how quickly you get there. Thus, 
the chaos reflected in the fractal weave of mixmaster uni- 
verses is unambiguous. 

It would be interesting to extend our study to include 
inhomogeneous collapse, and verify the connection be- 
tween temporal chaos and spatial turbulence (24| . As a fi- 
nal comment, we note that the chaos seen in the mixmas- 
ter system occurs at large curvatures. As is well known, 
most of the oscillations and bounces happen after Planck 
scale curvatures have been reached, so quantum effects 
cannot be ignored. Nonetheless, our results are consis- 
tent with the contention that generic classical singular- 
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